#ifndef POSITION_MODEL_H
#define POSITION_MODEL_H
#include "dur_hmm.h" 
#include <vector>
#include <armadillo>

class PositionModel
{
 public:
  typedef arma::vec dvec_t;
  typedef std::vector<dvec_t> vvec_t;
  typedef std::vector<int> sequence_t;
  typedef struct Variables {
    vvec_t eprob;
    dvec_t pprob;

    vvec_t ACC_eprob;
    dvec_t z, ACC_pprob; // z = prob(pos|seq)
    double seq_prob;
  } variable_t;
  typedef std::vector<sequence_t> vseq_t;
 public:
  int motif_length_;        // length of motif
  int max_seq_length_;
  int n_observation_;  // number of emission values
  variable_t var_;
 public:
  PositionModel(int L, int n_obs, int max_length);
  void rand(variable_t& var);
  void print_params() const;
  void train(const vseq_t& data, int max_iter = 15);
  int map_pos(const sequence_t& seq);
 protected:
  void start_accumulate(variable_t& var);
  void compute_z(const sequence_t& seq, variable_t& var);
  void accumulate(const sequence_t& seq, variable_t& var);
  void reestimate(variable_t& var);
};

#endif
